Adaptive subtraction of multiples using the L1-norm

A strategy for multiple removal consists of estimating a model of the multiples and then adaptively subtracting this model from the data by estimating shaping filters. A possible and efficient way of computing these filters is by minimizing the difference or misfit between the input data and the filtered multiples in a least-squares sense. Therefore, the signal is assumed to have minimum energy and to be orthogonal to the noise. Some problems arise when these conditions are not met. For instance, for strong primaries with weak multiples, we might fit the multiple model to the signal (primaries) and not to the noise (multiples). Consequently, when the signal does not exhibit minimum energy, we propose using the L1-norm, as opposed to the L2-norm, for the filter estimation step. This choice comes from the well-known fact that the L1-norm is robust to ‘large’ amplitude differences when measuring data misfit. The L1-norm is approximated by a hybrid L1/L2-norm minimized with an iteratively reweighted leastsquares (IRLS) method. The hybrid norm is obtained by applying a simple weight to the data residual. This technique is an excellent approximation to the L1-norm. We illustrate our method with synthetic and field data where internal multiples are attenuated. We show that the L1-norm leads to much improved attenuation of the multiples when the minimum energy assumption is violated. In particular, the multiple model is fitted to the multiples in the data only, while preserving the primaries. I N T R O D U C T I O N A classical approach to attenuating multiples consists of building a multiple model (see e.g. Verschuur, Berkhout and Wapenaar 1992; Berkhout and Verschuur 1997), and adaptively subtracting this model from the data, which is contaminated with multiples, by estimating shaping filters (Dragoset 1995; Liu, Sen and Stoffa 2000; Rickett, Guitton and Gratwick 2001). The estimation of the shaping filters is usually carried out in a least-squares sense making these filters relatively easy to compute. By using the L2-norm, we implicitly assume that the resulting signal, after the filter estimation step, is orthogonal to the noise and has minimum energy. These assumptions might not hold, and other methods, such as pattern-based approaches (Spitz 1999; Guitton et al. 2001), have been proposed to avoid these limitations. For instance, ∗E-mail: [email protected] when a strong primary is surrounded by weaker multiples, the multiple model will match both the noise (multiples) and the signal (primaries), such that the difference between the data and the filtered multiple model is a minimum in a least-squares sense. Consequently, some primary energy might leak into the estimated multiples and vice versa. We therefore need to find a new criterion or norm for the filter estimation step. We propose estimating the shaping filters with the L1norm instead of the L2-norm, thus removing the necessity for the signal to have minimum energy. This choice is driven by the simple fact that the L1-norm is robust to ‘outliers’ (Claerbout and Muir 1973) and large amplitude anomalies. Because the L1-norm is singular where any residual component vanishes, we use a hybrid L1/L2-norm that we minimize with an iteratively reweighted least-squares (IRLS) method. This method is known to give an excellent approximation of the L1-norm (Gersztenkorn, Bednar and Lines 1986; Scales and Gersztenkorn 1987; Bube and Langan 1997; Zhang, C © 2004 European Association of Geoscientists & Engineers 27 28 A. Guitton and D.J. Verschuur Chunduru and Jervis 2000). The main property of the hybrid norm is that it is continuous and differentiable everywhere, while being robust for large residuals. In the first section, we illustrate the limitations of the leastsquares criterion using a simple 1D problem. We then introduce our proposed approach, based on the L1-norm, to improving the multiple attenuation results. In a second synthetic example, we attenuate internal multiples with the L2and L1-norms. Finally, we apply shaping filters to a multiplecontaminated gather from a seismic survey, showing that the L1-norm leads to substantial attenuation of the multiples. P R I N C I P L E S O F L 1N O R M A N D L 2N O R M S U B T R A C T I O N In this section, we demonstrate with a 1D example that the attenuation of multiples with least-squares adaptive filtering is not effective when strong primaries are located in the neighbourhood. This simple example leads to a better understanding of the behaviour of our adaptive scheme in more complicated cases. Shaping filters and the L2-norm In Fig. 1, a simple 1D problem is considered. Figure 1(a) shows four events corresponding to one primary (on the left) and three multiples (on the right). Note that the primary has a larger amplitude than the multiples. Figure 1(b) shows a mulFigure 1 (a) The data with one primary at 0.06 s and three multiples at 0.14 s, 0.2 s and 0.32 s. (b) The multiple model that we want to adaptively subtract from (a). tiple model that corresponds exactly to the real multiples. For L2-norm subtraction the goal is to estimate a shaping filter f that minimizes the objective function, e2(f) = ‖d − Mf‖2, (1) where M is the matrix representing the convolution with the time series for the multiple model (Fig. 1b) and d is the time series for the data (Fig. 1a). If we estimate the filter f with enough degrees of freedom (enough coefficients) to minimize (1), we obtain the estimated primaries, i.e. d − Mf (Fig. 2a), and the estimated multiples, i.e. Mf (Fig. 2b). The estimated primary signal does not resemble the primary in Fig. 1(a). In Fig. 3, the corresponding shaping filter is shown. Note that this filter is not a unit spike at lag = 0 as expected. The problem stems from the least-squares criterion which yields an estimated signal that, by definition, has minimum energy. In this 1D case, the total energy in the estimated signal (Fig. 2a) is e2 = 2.4, which is less than the total energy of the primary alone (e2 = 4). This is the fundamental problem if we use the L2-norm to estimate the shaping filter. In the next section, we show that it is better to use the L1-norm if the multiples and the primaries are not orthogonal in the L2-norm sense. Shaping filters and the L1-norm The strong primary in Fig. 1 can be seen as an outlier that receives much attention during the L2 filter estimation. C © 2004 European Association of Geoscientists & Engineers, Geophysical Prospecting, 52, 27–38 Adaptive subtraction of multiples using the L1-norm 29 Figure 2 (a) The signal estimated with the L2-norm. (b) The noise estimated with the L2-norm. Figure 3 Shaping filter estimated for the 1D problem with the L2-norm. One lag is equivalent to one time sample (0.004 s). This filter is not a single spike at lag = 0. The maximum value of the filter is one at zero lag and the minimum value is −0.2 at lags −5 and +9. Consequently, some of the signal we want to preserve leaks into the noise. Because the L1-norm is robust to outliers, we propose using it to estimate the filter coefficients. This insensitivity to a large amount of ‘noise’ has a statistical interpretation: robust measures are related to long-tailed density functions in the same way that L2 is related to the short-tailed Gaussian density function (Tarantola 1987). In this section, we show that the L1-norm solves the problem referred to in the preceding section. Our goal now is to estimate a shaping filter f that minimizes the objective function, e1(f) = |d − Mf|1. (2) C © 2004 European Association of Geoscientists & Engineers, Geophysical Prospecting, 52, 27–38 30 A. Guitton and D.J. Verschuur Figure 4 (a) The signal estimated with the L1-norm. (b) The noise estimated with the L1-norm. Figure 5 Shaping filter estimated for the 1D problem with the L1-norm. One lag is equivalent to one time sample (0.004 s). This filter is a single spike at lag = 0 with amplitude value 1.0. The function in (2) is singular where any residual component vanishes, implying that the derivative of e1(f) is not continuous everywhere. Unfortunately, most of our optimization techniques, e.g. conjugate-gradient or Newton methods, assume that the first derivative of the objective function is continuous in order to find its minimum. Therefore, specific techniques have been developed either to minimize or to approximate the L1-norm. For instance, various approaches based on linear programming have been used with success (see e.g. Barrodale and Roberts 1980). Other robust measures, such as the Huber norm (Huber 1973), can also be considered with an appropriate minimization scheme (Guitton and Symes 1999). Alternatively, our implementation is based on the minimization of a hybrid L1/L2-norm with an iteratively reweighted least-squares (IRLS) method (Gersztenkorn et al. 1986; Scales C © 2004 European Association of Geoscientists & Engineers, Geophysical Prospecting, 52, 27–38 Adaptive subtraction of multiples using the L1-norm 31 Figure 6 (a) A synthetic shot gather containing many internal multiples. (b) The internal multiple model, which exactly matches the internal multiples in (a). Figure 7 Histograms of the input data (Fig. 6a) and of the noise (Fig. 6b). The multiples have a much weaker amplitude distribution and the L1-norm should be used. and Gersztenkorn 1988; Scales et al. 1988; Bube and Langan 1997). This technique is known to give a good approximation of the L1-norm. In this case, the objective function we minimize becomes e1(f) = ‖W(d − Mf)‖2, (3)